Leaflet Exploratory

Load packages

pacman::p_load('tidyverse', 
               'sf', 
               'here',
               'janitor',
               'tmap', 
               'kableExtra',
               'patchwork',
               'ggthemes',
               'colorRamps',
               'stars',
               'terra',
               'leaflet',
               'leaflet.extras')

Read Rast (calfire)

rast_path <- here::here('data','ds1327', 'tiff', 'ds1327.tif')

gdb_path <- here::here('data','ds1327', 'ds1327.gdb')

habitat <- rast(gdb_path)

veg_type <- habitat

activeCat(veg_type, layer = 1) <- 'LIFEFORM'

veg_type <- project(veg_type, "EPSG:4326") %>% 
  aggregate(fact = 20)

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

Read GAP GDB

gap <- st_read(here::here('data',
                          'PADUS4_1_State_CA_GDB_KMZ', 
                          'PADUS4_1_StateCA.gdb'), 
                          quiet = TRUE) %>% 
  clean_names() 

Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
automatically selected the first layer in a data source containing more than
one.

Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
Message 1: organizePolygons() received a polygon with more than 100 parts. The
processing may be really slow.  You can skip the processing by setting
METHOD=SKIP, or only make it analyze counter-clock wise parts by setting
METHOD=ONLY_CCW if you can assume that the outline of holes is counter-clock
wise defined

sf_use_s2(FALSE)

Spherical geometry (s2) switched off

gap_clean <- gap %>% 
  select(own_type, gap_sts, mang_type, Shape) %>%
  st_cast("MULTIPOLYGON") %>%  # Convert curves FIRST
  st_make_valid() %>%
  st_transform(4326) %>%
  st_simplify(dTolerance = 0.01) %>%
  group_by(gap_sts) %>%
  summarise(geometry = st_union(Shape), .groups = 'drop')

Warning in st_simplify.sfc(st_geometry(x), preserveTopology, dTolerance):
st_simplify does not correctly simplify longitude/latitude data, dTolerance
needs to be in decimal degrees

although coordinates are longitude/latitude, st_union assumes that they are
planar

although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar
although coordinates are longitude/latitude, st_union assumes that they are
planar

sf_use_s2(TRUE)

Spherical geometry (s2) switched on

Read Point Count

point_count <- read_csv(here::here('data', 'point_count.csv')) %>% 
  clean_names()

New names:
Rows: 683669 Columns: 39
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," chr
(19): GlobalUniqueIdentifier, ProjectCode, ProjectName, LocalityID, Stu... dbl
(16): ...1, SamplingUnitId, ParentSamplingUnitId, DecimalLatitude, Deci... lgl
(2): ProportionAreaSurveyed, InFocalHabitat date (1): ObservationDate time (1):
Time
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`

Read area search

area_search <- read_csv(here::here('data', 'area_search.csv')) %>% 
  clean_names()

Rows: 163889 Columns: 37
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (15): GlobalUniqueIdentifier, ProjectCode, ProjectName, LocalityID, Stu...
dbl  (15): SamplingUnitId, ParentSamplingUnitId, DecimalLatitude, DecimalLon...
lgl   (5): AreaSurveyed, DistanceFromObserver, Salinity, Temperature, DepthMin
date  (1): ObservationDate
time  (1): Time

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Join both and make Geodatabase

point_area_geo <- full_join(area_search, point_count) %>% 
  st_as_sf(coords = c("decimal_longitude", "decimal_latitude"), crs = 4326) 

Joining with `by = join_by(global_unique_identifier, project_code,
project_name, locality_id, study_area, sampling_unit_id,
parent_sampling_unit_id, sampling_unit_type_name, decimal_latitude,
decimal_longitude, visit, protocol_code, observation_date, year_collected,
month_collected, day_collected, julian_day, julian_day_ve, time, collector,
scientific_name, common_name, species_code, phylogen_order, taxon_id,
distance_from_observer, observation_count, no_observations,
record_permissions)`

Leaflet layers

# Create color palettes
gap_pal <- colorFactor(
  palette = c("green", "yellow", "orange", "red"),  # Adjust colors as needed
  domain = gap_clean$gap_sts  # Or whatever your column name is
)



# Map with colored polygons
leaflet_map <- leaflet() %>% 
  addProviderTiles(providers$Esri.WorldTerrain, group = "ESRI Terrain") %>% 
  addMiniMap(toggleDisplay = TRUE) %>% 
  
  # GAP Status with colors
  addPolygons(data = gap_clean, 
              group = 'GAP Status',
              fillColor = ~gap_pal(gap_sts),  # Replace gap_sts with your column
              color = "black", 
              weight = 1, 
              fillOpacity = 0.6) %>% 
  
  # Vegetation Type with colors
  addRasterImage(veg_type,
                 group = 'Vegetation Type',
                 opacity = 0.8) %>% 
  
  # Add legends
  addLegend(pal = gap_pal, 
            values = gap_clean$gap_sts,  # Replace with your column
            title = "GAP Status",
            position = "bottomleft",
            group = 'GAP Status') %>%
  
  addLayersControl(
    overlayGroups = c('GAP Status', 'Vegetation Type'),
    options = layersControlOptions(collapsed = TRUE)
  ) %>% 
  
  leaflet.extras::addResetMapButton()

leaflet_map